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The interplay between the structural relaxation and the rheological response of a simple amorphous system 
(an 80:20 binary Lennard- Jones mixture [W. Kob and H.C. Andersen, PRL 73, 1376 (1994)] is studied via 
molecular dynamics simulations. In the quiescent state, the model is well known for its sluggish dynamics 
and a two step relaxation of correlation functions at low temperatures. An ideal glass transition temperature of 
T c = 0.435 has been identified in the previous studies via the analysis of the system's dynamics in the frame 
work of the mode coupling theory of the glass transition [W. Kob and H.C. Andersen, PRE 51, 4626 (1995)]. 
In the present work, we focus on the question whether a signature of this ideal glass transition can also be 
found in the case where the system's dynamics is driven by a shear motion. Indeed, the following distinction in 
the structural relaxation is found: In the supercooled state, the structural relaxation is dominated by the shear 
at relatively high shear rates, j, whereas at sufficiently low 7 the (shear-independent) equilibrium relaxation 
is recovered. In contrast to this, the structural relaxation of a glass is always driven by shear. This distinct 
behavior of the correlation functions is also reflected in the rheological response. In the supercooled state, the 
shear viscosity, 77, decreases with increasing shear rate (shear thinning) at high shear rates, but then converges 
toward a constant as the 7 is decreased below a (temperature-dependent) threshold value. Below T c , on the 
other hand, the shear viscosity grows as r\ cx I/7 suggesting a divergence at 7 = 0. Thus, within the accessible 
observation time window, a transition toward a non-ergodic state seems to occur in the driven glass as the driving 
force approaches zero. As to the flow curves (stress versus shear rate), a plateau forms at low shear rates in the 
glassy phase. A consequence of this stress plateau for Poiseuille-type flows is demonstrated. 

PACS numbers: 64.70.Pf,05.70.Ln,83.60.Df,83.60.Fg 



I. INTRODUCTION 

Although simple at first sight, suspensions of 
spherical colloidal particles under shear exhibit a 
rich phenomenology. In the dilute regime, at tem- 
peratures corresponding to the liquid state, forced 
Rayleigh scattering experiments 1 show an increase 
of diffusion constant upon shearing (shear thin- 
ning), distinct from Taylor dispersion 2 . Indeed, 
Taylor dispersion (displacement of particles along 
the flow direction, they move in the direction 
of shear gradient) would give rise to an enhance- 
ment of mobility via (Ax 2 (t)) ~ t 3 , whereas shear 
thinning leads to (Ax 2 (t)} = D{p/)t, with a dif- 
fusion coefficient, D(^), which increases with in- 
creasing 7 in a non-trivial fashion. 

When shearing starts in the crystalline phase, 
however, shear induced melting of the crys- 
talline structure is observed both in light scatter- 
ing experiments^* 4 - as well as in Brownian dynam- 
ics simulations- Furthermore, in the shear melted 
regime, experiments show evidence for shear thin- 
ning due to the presence of freely slipping two di- 
mensional crystalline layers^ 4 -. Simulations also 
show a shear thinning regime below a "critical" 
shear rate, j c , followed by a transition to a string- 
like order for 7 > 7^. 



On the other hand, studies of disordered suspen- 
sions of hard spheres show that shear thinning and 
shear melting phenomena may also occur in the ab- 
sence of a crystalline structur o 6 ' 7 ' 8 ' 9 . Similar obser- 
vations have also been made in light scattering echo 
studies of (disordered) dense emulsions^. Brown- 
ian dynamics simulations show that shear thinning 
in concentrated colloidal suspensions is related to 
the fact that, in the limit of low shear rates, the main 
contribution to the shear stress originates from the 
Brownian motion of colloidal particles and that this 
contribution decreases with 7 11 . 

Many theoretical approaches have been proposed 
in recent years aiming at a description of the rhe- 
ology of disordered materials. Some examples are 
studies of driven p-spin glasses^, the so called 
"soft glassy rheology" model (SGR) 13 and non- 
equilibrium extensions of the mode coupling theory 
(MCT) of the glass transitio n 14 ' 15 ' 16 ' 17 . 

These approaches have the common feature of 
predicting (or, more precisely, "reproducing", as the 
phenomenon was known before the theories have 
been developed) the shear thinning phenomenon in 
an ergodic disordered system, i.e. in a system where 
temporal correlations decay to zero in a finite time. 
In these systems, shear thinning is commonly at- 
tributed to the competition between a typical struc- 
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tural relaxation time (of the quiescent system) and 
the time scale imposed by the external drive. A 
shear thinning behavior is expected as soon as the 
latter becomes significantly shorter than the former, 
in which case the structural relaxation is said to be 
driven by the external force 12 . This idea is also sup- 
ported by the observation that even a simple liquid 
may exhibit shear thinning provided that the shear 
rate is high enough in order for shear disturbance to 
dominate the system dynamics^. 

However, the statements of different theories di- 
versify as it comes to the rheological response of a 
"jammed system" such as an amorphous solid. The 
discrepancy is best seen when considering the be- 
havior of the yield stress, i.e. the stress response 
in the limit of vanishing shear rate (driving force), 
cr yield = <r(7 — ¥ 0). While calculations based 
on the p-spin glasses in the thermodynamic limit 
predict that no yield stress exists in the system 12 
(c"yieid vanishes with a power law, cr yie i d oc 7" with 
< a < 1), the SGR model 13 predicts the onset 
of a (dynamic) yield stress at the "jamming transi- 
tion", (Tyieid = 1 — x. Here, x is a noise temper- 
ature which controls the "degree of jamming" or 
the distance from the glass transition, x = 1 corre- 
sponds to the glass transition (or "jamming") tem- 
perature, and < x < 1 characterizes the glassy 
or "jammed" phase (for noise temperatures above 
the transition, 1 < x < 2, a power law decrease 
of the shear viscosity with the applied shear rate is 
found-^). 

Interestingly, the presence of a finite yield stress 
in an amorphous solid has also been predicted 
within the simplest (or idealized) version of the 
non-equilibrium MCT approach 14,1 -. A related 
MCT approach to the fluctuations around the steady 
state has recently been proposed in reference^^. 
The issue of yield stress, however, could not be ad- 
dressed in that approach. 

Despite the above mentioned qualitatively differ- 
ent predictions on the stress response of a glass, 
both the p-spin approach and the non-equilibrium 
MCT commonly predict that the presence of an ex- 
ternal drive leads to a melting of the amorphous 
solid. 

In reference 12 it is argued that the main reason for 
the absence of a yield stress in the p-spin model lies 
in its mean-field nature so that free energy barriers 



are impenetrable. If the system is quenched from 
a high temperature to a temperature in the glassy 
phase while at the same time being driven by a non- 
conservative force, it remains sliding above the free 
energy threshold below which the free energy sur- 
face is split into exponentially many disconnected 
regions 12 . This threshold does not exist at high 
temperatures but emerges at a critical temperature, 
T c , separating the liquid like region from the glassy 
phase. On the other hand, if the system is prepared 
in an energy state above the free energy threshold, it 
starts evolving towards the threshold free energy as 
the temperature is reduced below T c . This process 
slows down continuously due to the decreasing con- 
nectivity of the visited free energy landscape so that 
the free energy threshold is never reached. There- 
fore, even a vanishingly small driving force is suf- 
ficient in order to stop aging and keep the system 
sliding over this threshold free energy^ 2 -. 

This picture motivated a later work, where the ex- 
istence of a yield stress could be shown for a p-spin 
system with a finite number of spins 19 . The work 
follows the idea that free energy barriers are finite at 
finite system size, thus allowing the thermal activa- 
tions to play a role which is not possible in the case 
of an infinite system. In a system with finite free 
energy barriers, thermal activation becomes impor- 
tant since it allows the system to leave local free 
energy minima where it is trapped so that deeper 
valleys can be reached. A finite force must then 
be applied to the system in order to stop this aging 
process and to impose a steady sliding motion. For 
p = 3, the finite-sized version of the model has been 
investigated by Monte Carlo simulations supporting 
the existence of a critical driving force below which 
the system is trapped ('solid') and above which it 
"flows" ('liquid')^. 

A manifest of shear melting of an amorphous 
solid can be seen in the behavior of particle dis- 
placements. While in the quiescent state a particle 
is practically eternally trapped in its local environ- 
ment (the 'cage" formed by its nearest neighbors) 
it can explore far larger distances as the system is 
exposed to an external drive such as a shear motion. 
Thus, the particles in a driven glass are not localized 
but behave as in a fluid where the whole available 
space can be explored via diffusive motion. 

This "freedom" in particle motion is also reflected 
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in the behavior of the structural correlations. While 
in the quiescent system the structural relaxation 
times continuously grow as the system is quenched 
into a glassy phase (aging), eventually exceeding 
available experimental time window, all correla- 
tions decay to zero within a finite time in a driven 
glass. 

Based on the above picture of shear-driven melt- 
ing of an amorphous solid, one could possibly ex- 
pect that, regardless of the shear rate, a finite stress 
must always be applied in order to enable the sys- 
tem to overcome free energy barriers so that struc- 
tural relaxation can take place. If this conclusion is 
true, a stress plateau is expected as soon as the time 
scale of the inherent structural relaxation is far be- 
yond the accessible simulation time, since, in this 
case, the finite relaxation time in the driven system 
can be considered as uniquely imposed by the exter- 
nal drive, a situation very similar to shear-induced 
melting. 

In order to test this idea, we performed large scale 
molecular dynamics simulations of a simple molec- 
ular glass (see below for a description of the model) 
focusing on the interplay between the structural re- 
laxation and the rheological response of the sys- 
tem both in the supercooled state and in the glassy 
phase. As published in a recent paper—, our sim- 
ulation results clearly support the existence of a 
stress plateau at low shear rates in an amorphous 
solid. The focus of that paper was on the stress 
response and on a comparison of simulated results 
with MCT-based theoretical calculations. Here, we 
go a step further and study the dynamics of struc- 
tural relaxation in the same range of temperatures 
and shear rates for which the flow-curves have been 
studied in reference 20 . Furthermore, we will inves- 
tigate an interesting consequence of the presence of 
a stress plateau (which in this case plays the role of 
a yield stress) on a Poiseuille-type flow in a planar 
channel. It will be shown that in a region around the 
channel center where the stress is below the plateau 
stress the system behaves as an amorphous solid 
whereas it exhibits liquid-like behavior beyond this 
central region. 

It is noted that the rheological response of ex- 
actly the same model has already been studied in 
reference 2 ^-. However, due to a rather restricted 
range of studied shear rates, no conclusion could 



be made on the existence or not of a dynamic yield 
stress in the glassy phase. The novelty of our simu- 
lations is the extended range of shear rates allowing 
a conclusive statement on the existence of a stress 
plateau for the present model in the glassy state. 
When restricted to the range of shear rates acces- 
sible to our simulations, this stress plateau may be 
regarded as the yield stress of the system. 

The paper is organized as follows. After an in- 
troduction of the model and the simulation method 
in the next section, the effect of shear on the struc- 
tural relaxation of an amorphous system will be dis- 
cussed in section |nll In particular, it will be shown 
that the structural relaxation of an amorphous sys- 
tem exposed to a shear exhibits signature of the qui- 
escent state so that a study of the structural relax- 
ation alone is sufficient in order to find out whether 
the quiescent state of the (driven) system belongs 
to a liquid-like (supercooled) state or to a glassy 
phase. In section [IV] the rheological response of 
the system is studied. The focus of this section is 
to show that a similar information on the quiescent 
state of the system may also be gained via a study 
of the stress response. This is possible since the 
stress-shear rate curves of a glass are qualitatively 
different compared to those of a supercooled liquid. 
The information on the glass transition is, therefore, 
also encoded in the rheological response. Section IVl 
compiles our results. 



H. MODEL AND SIMULATION METHOD 

A generic glass forming system, consisting of 
an 80:20 binary mixture of Lennard- Jones particles 
(whose types we call A and B) at a total density of 
P = Pa+Pb = 1-2 and in a cubic box of length L = 10 
(N = 1200 particles) is studied 2 ^. 

A and B particles interact via t/ LJ (r) = 
4e af 3[(d af3 /r) 12 - (d af3 /r)% with a,/3 = A,B, e AB = 
1-5caa> 6 bb = 0.5e AA , c/ AB = 0.8<iAA> <^bb = 0.88<iAA 
and tub = m A . The potential was truncated at 
twice the minimum position of the LJ potential, 
r c ,a/3 = 2.2Abd a fj. The parameters e A A, <^aa and 
m A define the units of energy, length and mass. All 
other quantities reported in this paper are expressed 
as a combination of these units. The unit of time, 
for example, is given by t L j = d AA ^m A /e AA and 
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that of stress by €aa/ *^aa - Equations of motion are 
integrated using a discrete time step of dt = 0.005. 

The system density is kept constant at the value 
of 1.2 for all simulations whose results are reported 
here. This density is high enough so that no voids 
occur at low temperatures and low enough so that 
system dynamics remains sensitive to a variation of 
temperature (see references" for effects of high 
density/pressure on the liquid-glass transition). 

The present model was found suitable for an anal- 
ysis of many aspects of the mode coupling theory 
of the glass transition". In particular, at a to- 
tal density of p = 1.2, equilibrium studies of the 
model showed that the growth of the structural re- 
laxation times at low temperatures could be ap- 
proximately described by a power law as predicted 
by the ideal MCT, r relax oc (T — T C )~ 7MCT . Here, 
T c = 0.435 is the mode coupling critical temperature 
of the model and 7mct is the critical exponent. For 
the present binary Lennard- Jones system, numeri- 
cal solution of ideal MCT equations yields a value 
of 7mct ~ 2.5 27 . A similar value is also obtained 
for a binary mixture of soft spheres 28 . 

Simulation results are averaged over 10 inde- 
pendent runs. For this purpose, ten independent 
samples are equilibrated at a temperature of T = 

0. 45 (above T c ) and serve as starting configura- 
tions for all simulated temperatures and shear rates. 
The temperature is controlled via Nose-Hoover 
thermostat". It is set from T = 0.45 to the de- 
sired value at the beginning of shear, whereby only 
the ^-component of particle velocities is coupled to 
the heat bath (x being the streaming and z the shear 
gradient directions, see also below). 

The temperature quench is done only in one step, 

1. e. without a continuous variation from T start to 
T cnd . However, as the numerical value of T is 
changed, it takes a time of the order of the velocity 
autocorrelation time for the new temperature to be 
established. During this period of time the Maxwell 
distribution of velocities undergoes changes in or- 
der to adapt itself to a distribution determined by 
the new temperature. This time is of order unity (in 
reduced units) and quite short compared to all other 
relevant timescales in the problem. 

Previous studies of the stress-strain relation of the 
same model showed that the initial transient behav- 
ior is limited to strains below 50%2i.. Indeed, by 



shifting the time origin in measurements of vari- 
ous correlation functions, we verified that the time 
translation invariance was well satisfied in sheared 
systems for strains larger than 50%. We neglected 
strains 7 < 100% before starting the measurements. 
Unless otherwise stated, all simulations reported 
below had a length corresponding to a strain of 7.8 
(780%). In the steady state, correlation functions 
were averaged both over independent runs and over 
time origins distributed equidistantly along each 
simulation run. The shear stress is calculated using 
the virial expression^ 

-ziAr if p \ 

@ — ®xz ( / j T^i^ix^iz 2 / j •Eij-T zij ) > 

\ i i+j I 

(1) 

where (...) stands for statistical averaging, Vj and 
mi are the velocity and the mass of i-th particle, 
Xij = Xj — Xi and F zi j the z-component of the force 
of particle j on 2. 

Recent studies of the present model in the glassy 
state showed that the system may exhibit shear- 
localization if the shear rate is imposed by using 
a conventional Couette cell with moving atomistic 
walls'. A shear banding is, however, undesired in 
the context of present analysis since we are inter- 
ested in the effects of a homogeneous shear. On 
the other hand, simulations of the present model 
using the so-called Lees-Edwards boundary condi- 
tions along with the SLLOD equations of motion 
first proposed by Evans and coworkers 1 - show that a 
linear velocity profile forms across the system thus 
leading to a spatially constant velocity gradient'. 
We, therefore, also used this approach for our sim- 
ulations (see Fig. Q]for an illustration of the Lees- 
Edwards boundary condition). Within this simula- 
tion method, we do indeed observe a spatially con- 
stant velocity gradient in all studied cases (Fig. [2]). 

The external shear does work on the system. 
Therefore, in order to keep the system temperature 
at a prescribed value, this extra heat must be re- 
moved. This is done by coupling the y-component 
of particle velocities to the so-called Nose-Hoover 
thermostat''. We impose a flow in the x-direction 
with a shear gradient in the ^-direction. The y- 
component of particle velocities are therefore not 
affected by the local flow velocity, thus simplify- 
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vt dy/dt=V/L 



to the following set of equations of motion, 
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FIG. 1 : An illustration of the Lees-Edwards boundary condition. The 
central box depicts a simulation box containing two particles only. 
All surrounding cells are images of the central box. As indicated by 
bold arrows, the image box on the top is moved with a constant ve- 
locity of V towards right, whereas the image box on the bottom is 
moved with the same velocity in the opposite direction. This defines 
an overall shear rate of j — V/L, where L is the box length in the 
vertical direction. Note that all particles of a given image box have 
exactly the same additional velocity V (—V). The method is there- 
fore sometimes referred to as the sliding brick boundary condition. 
This image is motivated by a similar figure in reference 18 . 
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FIG. 2: Velocity profiles across the system obtained for various tem- 
peratures as indicated (ranging from the glassy phase to the super- 
cooled regime). As seen here, shear localization is absent at all stud- 
ied temperatures. However, at the global shear rate studied here, 
shear bands may occur if the shear is imposed via the motion of atom- 
istic walls^. The use of the SLLOD algorithm in conjunction with 
the Lees-Edwards boundary condition is, therefore, essential in order 
to ensure a homogeneous shear across the system. 



ing the computation of instantaneous kinetic energy 
which enters the thermostating part of equations of 
motion. This choice also ensures that possible arte- 
fact due to a profile biased thermostat 33 are absent. 
Further details about the simulation and the model 
can be found in reference 31 . 
The above choice of the simulation method leads 
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Here {x u y h z*), (p xi , p yi , p zi ) and denote 
coordinates, momenta and the mass of i-th particle 
respectively. N = N A + N B is the total number of 
particles, k B the Boltzmann constant and F the total 
force on a particle. It is noted that the momenta oc- 
curring in the above equations are the peculiar ones, 
i.e., they correspond to the particles' momenta in 
a (local) frame of reference moving with the flow 
(ip) — 0) 18 . The variable, £, plays the role of a 
friction (acceleration) coefficient, since a positive 
(negative) £ tends to decrease (increase) p y . The 
variation of £, on the other hand, is controlled by the 
deviation of the actual kinetic energy of the system 
from that prescribed by the temperature, T, of the 
heat bath. The parameter Q controls the strength of 
the coupling of the particles' momenta to the heat 
bath: the smaller Q the stronger the coupling (see 
reference 34 for more details on thermostating meth- 
ods). 

III. STRUCTURAL RELAXATION IN A DRIVEN SYSTEM 

It is one of the aims of the present section to show 
that the structural relaxation of an amorphous sys- 
tem exposed to a shear already exhibits signature of 
the quiescent state of the system. In other words, a 
study of the structural relaxation in the driven case 
allows one to find out whether the quiescent state 
of the system belongs to a liquid-like (supercooled) 
state or to a glassy phase. 

As will be shown below, this possibility of dis- 
tinguishing between a sheared glass and a driven 
supercooled liquid is closely related to the fact that, 
in the supercooled state, the structural relaxation is 
affected by the external drive at high shear rates 
only. As the shear rate tends to zero, the (shear- 
independent) equilibrium relaxation is recovered. 
In contrast to this, the structural relaxation of a 
sheared glass is always driven by the imposed shear, 
since, in an amorphous solid, there is no equilib- 
rium relaxation in the available experimental (or 
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simulation) time window. 

In order to illustrate these features, we show in 
Figs. [3] and |4] the incoherent scattering function (q y 
is the y-component of the wave vector), 

1 N 

$ «(0 = ^(E ex P^W - W(0))l), ( 2 ) 

i=l 

obtained from molecular dynamics simulations at 
various shear rates and temperatures. The use of 
y-component naturally eliminates undesired contri- 
butions to particle rearrangements arising, for ex- 
ample, from the motion of particles with the flow 
and the so-called Taylor dispersion (affecting the 
displacement of particles along the flow direction 
as they move in the direction of shear gradient) 2 -. 

We measure $ g (t) for a wave vector of q — 7.1, 
corresponding to the principal peak of the static 
structure factor. thus reflects collective par- 

ticle rearrangements on the scale of average near- 
est neighbor distance, the relevant length scale for a 
study of the cage effect. 

In Fig. [3l & g (t) is depicted for two temperatures 
from the supercooled regime. In each panel, the 
incoherent scattering function is shown for various 
shear rates ranging from a regime, where the sys- 
tem dynamics is accelerated by the imposed shear, 
to lower 7 where the structural relaxation is inde- 
pendent of shear rate. First note that, at short times, 
all curves overlap regardless of the imposed shear 
rate indicating that the short time dynamics is un- 
affected by the shear. The situation is, however, 
different when long time behavior of is con- 
sidered. At high 7, the final relaxation of $ ? (t) 
strongly depends on the magnitude of the imposed 
shear rate: The higher 7 the faster the decay of 
$ 9 (t). However, as 7 is progressively decreased, 
the dependence of $ g (t) on 7 becomes weaker, 
eventually disappearing for the lowest shear rates 
shown in the figure. 

A comparison of <& 9 (i) for a shear rate of 7 = 1CT 4 
in both panels of Fig. [3] shows that a shear rate 
which does not affect the equilibrium relaxation at 
a high temperature (here, T = 0.525) may induce 
an acceleration of the system dynamics at a lower 
temperature (T = 0.45). Interestingly, as will be 
shown in the next section (Fig. [8]), the stress re- 
sponse exhibits a linear (non linear) behavior for 



roughly the same shear rates where $ 9 (t) is unaf- 
fected (affected) by shear. 

The above comparison follows a similar plot for 
T = 0.5 in reference 2 ^ (see the upper left panel in 
Fig. 8). The new aspect in our plots is the study of 
considerably wider range of shear rates allowing to 
decrease 7 sufficiently in order to recover equilib- 
rium behavior at a temperature as low as T = 0.45, 
a feature suggested but not shown in reference 2 ^-. 

The discussion of Fig. [3] underlines the idea that 
the linear or non linear nature of the system re- 
sponse to an imposed shear rate is the outcome of 
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FIG. 3: Effect of a homogeneous shear rate on the relaxation be- 
havior of the incoherent scattering function. <fr q (t) is shown for two 
temperatures, both belonging to the supercooled regime. The overall 
behavior is similar at both temperatures: At high shear rates an ac- 
celeration of the system dynamics is observed, whereas equilibrium 
curves are recovered at sufficiently low 7. However, a comparison 
of these two panels clearly shows that the linear response regime is 
shifted towards considerably lower shear rates as temperature is re- 
duced. The mode coupling critical temperature of the system is T c = 
0.435^). From left to right 7 = 1CT 1 , 3 x 1(T 2 , 10~ 2 ,...,3x 
1(T 5 forT = 0.525 and 7 = Kr\ 3 x 1(T 2 , 1(T 2 , 3 x 1CT 6 
for T = 0.45. 
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the competition between the time scale imposed by 
the flow and the inherent (structural) relaxation time 
of the system. The faster the inherent system dy- 
namics, the larger the range of shear rates in the 
linear regim o 12ilSi16 . 

Figure @] illustrates the effect of shear on the re- 
laxation behavior of a glass. For this purpose, $ g (t) 
is depicted for a temperature far below the mode 
coupling critical temperature of the model. Again, 
at short times, all curves are identical irrespective of 
the imposed shear rate. As for the long time behav- 
ior, ® q (t) exhibits a dependence on 7 at all simu- 
lated shear rates. The plot in Fig.|4]is also motivated 
by the upper right panel in Fig. 8 of reference 21 . 
In addition to a wider range of shear rates shown 
in our plot, Fig. @] also compares the sheared 
to its quiescent counterpart, t w ; 7 = 0), mea- 
sured after a waiting time of t w = 10 5 (and aver- 
aged over independent initial samples). This allows 
to unambiguously demonstrate the fact that, even at 
the lowest studied shear rate (7 = 3 x 1CT 6 ), the 
final decay of the relaxation function at T = 0.3 is 
indeed a result of imposed shear. 

In the linear response regime, transport coeffi- 
cients can often be expressed as time integrals of 
equilibrium correlation functions. The shear vis- 
cosity, for example, is given by the well-known 
Green-Kubo relation rj = f °° g(t)dt, where g(t) = 
V/ {k^T)(a X y{t)a X y(Q)) (V denotes the system vol- 
ume and a xy a component of the stress tensor with 
zero mean, (<J xy ) = 0). Thus, in a range where g(t) 
becomes independent of the applied shear rate, the 
shear viscosity must also be a constant (Newtonian 
behavior). 

As can be seen from Eq. (Q~|), the shear stress 
can be expressed as a function of particle positions 
and (much faster variables) momenta and forces, 
whereas & q (t) depends on the particle positions 
only. Thus, at least for wave vectors corresponding 
to the average interparticle distance, one expects a 
faster decay of stress autocorrelation function com- 
pared to a relaxation of & q (t). This expectation is 
born out in Fig. [5J However, despite the fact that the 
time scale for the final decay of g(t) is roughly by 
one order of magnitude shorter than the decay time 
of (for q being the maximum of the struc- 

ture factor), the cross over from the shear thinning 
behavior to equilibrium relaxation occurs at practi- 




time [LJ unit] 

FIG. 4: Same as in Fig.[3]but for a temperature below T c . The inco- 
herent scattering function, & q (t), shows a pronounced two step re- 
laxation. Whereas the short time dynamics is practically unaffected 
by shear, the long time behavior of $ g (t) is clearly dominated by 
the imposed shear rate. The lower 7 the larger r, the time necessary 
for a substantial decay of & q (t). For comparison, the quiescent case, 
$ 9 (t , t„ ; 7 = 0) , measured after a waiting time of t w = 10 , is also 
shown (see diamonds; in the case of $ q (t,t w ;7 = 0); the averag- 
ing is done over a large number of initial configurations). Obviously, 
after a waiting time of i w = 10 , the inherent system dynamics is 
far slower than the shear driven decay of for all studied shear 
rates. This must be contrasted to the supercooled state, where equi- 
librium behavior is recovered at low 7 (Fig. [3). From left to right, 
7 = 10~\ 3 x 10~ 2 , 10~ 2 ,...,3 x 10" 6 . 

cally the same 7 for the both types of correlation 
functions (compare the data for T = 0.45 in Figs. [3] 
and©. 

This observation underlines the relevance of 
$> q (t) as an appropriate quantity for at least quali- 
tative studies of non-linear rheology. In a 7-range 
where & q (t) becomes independent of the shear rate, 
the shear viscosity is also expected to be indepen- 
dent of 7. The study of the flow curves presented 
below (Fig. [8]) confirms this expectation. 

We examine the validity of the so-called "time- 
shear superposition principle", predicted both 
within the spin glass theory^ and the non- 
equilibrium MCT— . This property indicates that the 
shape of & q (t, 7) at sufficiently large times is inde- 
pendent of shear rate. More precisely, $ g (t,7) = 
F q {t), where i = t/r^j) is a dimensionless time. 
The shear rate dependent time, t{^), characterizes 
the time scale of the final decay of the correlator 
(the so-called a-relaxation). 

This property is checked for in Fig. [6] In this fig- 
ure, the same data as shown in Fig. [3] are depicted 
versus rescaled time t/r(j) (r is estimated from 
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FIG. 5: Effect of a homogeneous shear rate on the relaxation be- 
havior of the stress autocorrelation function, (a xy (t)a xy (0)) for a 
temperature of T = 0.45 (supercooled state) and T = 0.3 (glassy 
phase). In the supercooled state, an acceleration of the stress relax- 
ation is observed at high shear rates, whereas equilibrium curves are 
recovered at sufficiently low j. In the glass, on the other hand, no 
equilibrium curve exists within the simulation time window so that 
even the smallest shear rate studied significantly alters the system dy- 
namics. From left to right 7 = 10" \ 3xl0~ 2 , 10~ 2 , 3 x 10" 6 . 
Note that, contrary to o xz whose average gives the shear stress in 
the system, {a xy ) = 0. This explains why the studied correlation 
function decays to zero and not to a finite positive value. 



®q{t = t) = 0.1). As seen from the both panels of 
Fig. |6l there is a group of curves following a mas- 
ter curve. These curves belong to low 7, whereas 
at higher 7 the time-shear superposition principle is 
violated. This is reminiscent of similar deviations 
observed in the non-driven supercooled state where 
the role of 7 is played by temperatur e 26 ' 35 . 

For the both temperatures shown in Fig. [61 there 
is an intermediate time window ("intermediate" in 
the sense that it is long compared to the initial de- 
cay time but short compared to the final decay time 
of $ q (t)), known as MCT-/3— relaxation regime, for 
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FIG. 6: A test of time-shear superposition principle. Same data as 
in Fig.[5]are plotted versus rescaled time t/rip/). At low 7 and for 
long times, curves follow a master curve whereas deviations from the 
master curve are observed for higher shear rates. The thick solid line 
is a fit using the von Schweidler law,Eq. QJ, with / 9 = 0.79, 6 = 0.49 
and A = 0.783 for T = 0.525 and / 9 = 0.79, 6 = 0.53 and A = 0.922 
for T = 0.45. Note that the von Schweidler law is predicted within 
equilibrium MCT as an approximation to the late /3-relaxation (see 
Refsi 2 ^ 2 ^ for a more detailed test of ideal MCT predictions of the von 
Schweidler law). 



which the master curve can be well fitted to the so- 
called von Schweidler law, 



%(t) = f q - A(t/r) 



(3) 



Here, the non ergodicity parameter, f q , determines 
the height of the plateau (to which the short time 
decay leads) and b is the so-called von Schweidler 
exponent. We find practically the same numerical 
value for f q and very close values for b at both 
temperatures. This is consistent with ideal MCT, 
which yields a theoretical justification of this em- 
pirical law 36 . Indeed, here we use the results of 
equilibrium MCT, since at low shear rates (where 
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the theory is supposed to apply) recovers its 
equilibrium behavior. 

For the glass, a similar plot as in Fig. [6] is shown 
in Fig. [7] Obviously, also in the glass, the low shear 
rate data confirm the validity of the time- shear su- 
perposition principle. It is worth mentioning that, 
Eq. © remains valid also in the case of a driven 
amorphous solid. However, in contrast to the case 
of the quiescent system, where the von Schweidler 
exponent, b, may depend on the system properties 
at equilibrium, the non-equilibrium MCT predicts 
6=1 for the case of a driven glass 1415 . Using 6=1, 
a fit to Eq. © is also shown in Fig. [7J As also 
expected from the enhanced solid-like character of 
the system at T = 0.3 compared to T = 0.45 and 
T = 0.525, the value obtained for the non ergodic- 
ity parameter, f q = 0.84, is significantly larger than 
f q = 0.79, obtained in the supercooled state (Fig.©. 

IV. RHEOLOGICAL RESPONSE 

In this section, the rheological response of the 
system to a spatially constant (homogeneous) shear 
rate is studied. It will be shown that the stress re- 
sponse in the supercooled state is qualitatively dif- 
ferent from that in the glass. In the supercooled 
regime, the non linear response at high shear rates 
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FIG. 7: A test of time-shear superposition principle in the glass. 
Data in Fig. 5] are plotted versus rescaled time t/rft). At low 7 
and for long times, curves follow a master curve whereas deviations 
from the master curve are observed for higher shear rates. The thick 
solid line is a fit to Eq. $3^ with an exponent of b = 1—. This gives 
a non ergodicity parameter, f q — 0.84, which is larger than the one 
obtained in the supercooled state. 



is always followed by a linear relation between the 
shear stress and the shear rate (resulting in a con- 
stant viscosity) as the shear rate is sufficiently low- 
ered. However, as temperature is decreased, the lin- 
ear regime is shifted towards progressively lower 
shear rates and eventually disappears in the glassy 
phase. In other words, in the glass, any shear rate, 
as low as it might be, leads to a non linear response. 

In order to elucidate this aspect, we investigate 
the behavior of the flow curves (shear stress versus 
shear rate) across an ideal glass transition. 

Simulation results on the stress response to a time 
independent shear (flow-curves) are shown in the 
top panel of Fig. [8] for temperatures ranging from 
the supercooled state to the glassy phase. In the su- 
percooled state, the shear thinning behavior at high 
shear rates crosses over to the linear response as 
7^0. However, as the system is cooled toward 
the glassy phase, the linear response is reached at 
progressively lower shear rates. In the glass, on 
the other hand, the shear stress becomes indepen- 
dent of 7 at low shear rates. For many practical 
purposes (see e.g. Fig. [9]) this stress plateau plays 
the role of a yield stress: If the system is subject to 
stresses below this plateau stress, it behaves like a 
solid body (zero velocity gradient) while for higher 
stresses it exhibits a liquid-like character. Let us 
mention that a stress plateau is also observed in re- 
cent experiments on the rheology of dense colloidal 
dispersion a 37 ' 38 . 

The bottom panel of Fig. [8] depicts the same 
data as in the top panel in a different way. Here, 
the apparent shear viscosity, 77 app = cr/j, is plot- 
ted against 7. The presence of a linear response 
regime in the supercooled state is now illustrated as 
a constant viscosity line, whereas the apparent I/7- 
divergence of the viscosity in the glass (for 7^0) 
is also nicely born out. From this panel, one can 
easily extract, for temperatures in the supercooled 
state, the shear rate at which the crossover from lin- 
ear to non linear response takes place. This gives 
7 = 3 x 10~ 4 for T = 0.525 and 7 = 3 x 10~ 5 for 
T = 0.45. Interestingly, these values also mark the 
beginning of the shear rate dependence of 3> g (t) for 
the corresponding temperatures (Fig. [3]). 

Let us work out an interesting consequence of 
the presence of a yield stress on the flow behav- 
ior of the system. For this purpose, we simulate a 
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FIG. 8: Top: Simulated shear stress versus shear rate 
for various temperatures ranging from the glassy state 
to the supercooled regime (from top to bottom: T — 
0.01, 0.2, 0.3, 0.34, 0.36, 0.38, 0.4, 0.42, 0.43, 0.44, 
0.45, 0.47, 0.5, 0.525, 0.55, 0.6). In the supercooled state, 
the shear thinning behavior at high shear rates crosses over to the 
linear response as 7 — > 0. With decreasing temperature, the linear 
response is, however, reached at progressively lower 7. In the glass, 
on the other hand, the shear stress becomes independent of 7 at 
low shear rates (stress plateau). Figure adapted from reference 20 . 
Bottom: Apparent shear viscosity, rj app = 0-/7 for exactly the 
same temperatures as shown in the case of the shear stress. Not 
unexpectedly, for 7 — > 0, rj app reaches a constant value in the 
supercooled state, whereas it grows as I/7 in the glass. 



The above discussed Poiseuille-type simulation is 
performed in the following way. We first equilibrate 
a system of size L x x L y x L z = 30 x 30 x 86 
(containing 92880 particles) at a temperature of 
T = 0.45 and then quench it to a temperature of 
T = 0.2. The system is then aged for a waiting 
time of t w = 10 4 in order to suppress the fast aging 
process which occurs directly after the T-quench 
and allow the system to solidify. After this period 
of time, two solid walls are introduced by immo- 
bilizing all particles whose ^-coordinate satisfies 
| z | > 40 (note that z e [-43 43] in general; this 
gives rise to walls of three particle diameter thick- 
ness). A flow is then imposed by applying on each 
particle a constant force of F e = 0.025 (in LJ units). 
In the case of a liquid at equilibrium, this would 
give rise to a parabolic velocity profile. However, 
as Fig. [9] clearly demonstrates, the velocity profile 
rather exhibits a behavior expected for a yield stress 
fluid: In the central region defined by the condi- 
tion of a < Oyjeid ~ 0.5 (see the stress plateau at 
T = 0.2 in the upper panel of Fig. [8]) the velocity 
profile is flat with a zero gradient while it gradu- 
ally departs from this constant behavior (shear rate 
becomming nonzero) beyond this central part of the 
channel. An interesting consequence of this behav- 
ior is that, if for a given driving force (pressure 
gradient) the channel width is too small to allow 
the formation of stresses above the system's yield 
stress, no flow will occur and the whole system will 
behave like a solid body. 



V. SUMMARY 



Poiseuille type flow in a planar channel. The study 
of such a situation is interesting since the stress in 
a Poiseuille-type flow is zero in the channel cen- 
ter and increases linearly with the distance from 
it. If the fluid under consideration exhibits (at least 
within the time window accessible for the simula- 
tion) a finite yield stress, one may expect that the 
fluid portion in a certain region around the chan- 
nel center (defined by the condition that the stress 
in this part of the channel be lower than the yield 
stress) should behave like a solid body while it 
should flow like a liquid further away from this re- 
gion. 



We report results of large scale molecular dy- 
namics simulations on the structural relaxation and 
the rheological response of a well established glass 
forming model, an 80:20 binary Lennard-Jones 
mixture first introduced by Kob and Andersen in the 
context of the dynamics of supercooled liquids 2 ^. 

Previous studies of the present mode l 22i25i2e 
showed that the model was suitable for an analysis 
of many aspects of the so-called mode coupling the- 
ory of the glass transition (MCT)^ 4 ^. In particular, 
for a total density of p = 1 . 2 , an ideal computer glas s 
transition (in the sense that the relaxation times ap- 
proximately obey a power law divergence predicted 
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transverse coordinate, z 

FIG. 9: Flow through a three dimensional planar channel of the 
model fluid studied in this paper. The flow is generated by imposing 
a constant body force of F c = 0.025 on each particle. At a tempera- 
ture of T = 0.2 and a total density of p — 1.2, the model exhibits a 
stress plateau at low shear rates which, for the present purpose, plays 
the role of a yield stress, <7 y i e id ~ 0.5 (see Fig. [8j. One can easily 
show that the stress in the channel behaves as a — pF c z, where z 
is the transverse coordinate (z — being the channel center)^. Ob- 
viously, for \z\ < 16.67(« cr y ieid/(pi ;ie )), the stress in the channel 
is below the yield stress. One thus expects the system to behave as 
a solid body for \z\ < 16.67: Either it should be at rest or move 
with a constant velocity (zero velocity gradient). Indeed, an inspec- 
tion of the velocity profile, u(z), and its derivative (rescaled to fit 
into the figure) confirms this expectation (the region delimited by the 
two vertical dashed lines). For \z\ > 16.67, on the other hand, the 
local shear stress exceeds the yield stress of the model system thus 
leading to its liquid like behavior. The system flows with a shear rate 
which non-linearly increases upon increasing stress. Note that if in 
the present case the channel width becomes smaller than w 33.34 
(Lennard- Jones units), there will be no region with stresses above 
Oyieid- As a consequence, there will be no liquid-like region. The 
whole system will behave like a solid body sticking to the bound- 
aries of the channel (flow blockage). This expectation is confirmed 
by our simulations (not shown here). 

by ideal MCT) was observed at a mode coupling 
critical temperature of T c = 0.435. 

The present work is strongly motivated by recent 
theoretical progress on the field of the rheology of 
dense amorphous systems. While various theoreti- 
cal approache o 121 ' 3 ' 15 provide a coherent description 
of the rheological response of a supercooled liquid 
(prediction of shear thinning at high shear rates (or 
driving force), 7, followed by linear response as 7 
approaches zero), they make different predictions 
with regard to the stress response of an amorphous 
solid at low shear rates. 

In particular, the issue of dynamic yield stress, 
Cyieid* (shear stress in the limit of 7 — > 0) remains 
controversial. While approaches based on the nu- 
merical studies of disordered p-spin systems in the 



thermodynamic limit suggest that cr y i eld identi- 
cally vanishes in the glassy phase, both the semi- 
phenomenological "soft glassy rheology" model 13 
(SGR, an extension of the trap model 41 taking into 
account the presence of an external drive) as well as 
the idealized version of the non-equilibrium MCT 
proposed in reference 14 predict the existence of a 
finite dynamic yield stress in the glassy phase. 

Nevertheless, qualitatively similar statements on 
the relaxation behavior of the correlation functions 
in a driven disordered solid are made both within 
the p-spin model and in the non-equilibrium MCT. 
In particular, shear melting of a glass and the recov- 
ery of the time translation invariance (shear stops 
aging 4 !) j s commonly predicted by both these ap- 
proaches. 

An explanation for the absence of dynamic yield 
stress in the p-sp'm model is given in reference^ re- 
lating this property to the presence of infinite free 
energy barriers occurring in the glassy phase in the 
thermodynamic limit (N — > 00). This picture has 
been checked by studying via Monte Carlo simula- 
tions a p-spin model with a finite number of spins, 
whereby considering a case with finite free energy 
barriers. Indeed, a yield stress in the glass is ob- 
served within these studies^. 

The effect of shear on the relaxation behavior of 
the intermediate scattering function, $ g (t), as well 
as the stress autocorrelation function, g(t), is stud- 
ied both in the supercooled state and in the glass. 
In the supercooled state, high shear rates lead to an 
accelerated decay of the correlation function (shear 
thinning) compared to the equilibrium relaxation 
which is recovered at sufficiently low shear rates. 

Interestingly, despite the fact that g(t) decays 
much faster than (for q corresponding to the 
inverse of the nearest neighbor distance), the cross 
over shear rate from shear thinning to equilibrium 
relaxation is practically the same for the both types 
of the correlation functions $ 9 (£) and g(t) (com- 
pare the data for T = 0.45 in Figs. [3] and [5]). Fur- 
thermore, this cross over shear rate also marks the 
change from linear to non-linear behavior in the 
stress response (Fig. [8]). These observations under- 
line the role of the intermediate scattering function 
as an appropriate observable for a study of at least 
qualitative features of non-linear rheology. 

In the glass, on the other hand, our simulations 
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clearly indicate that the dynamic yielding (shear 
melting) of an amorphous solid (Fig. [7]) is accom- 
panied by the presence of a finite stress plateau 
(Fig. [8]), which we identify as the yield stress of the 
model given the restriction to the accessible obser- 
vation time window. It is worth mentioning that a 
stress plateau is also found in recent experiments on 
the rheology of dense colloidal dispersion o 37 ' 38 . 

An interesting consequence of the presence of a 
yield stress is discussed for the case of a flow driven 
by an external body force (equivalent of a pressure 
driven Poiseuille-type channel flow; see Fig. HI). 
This choice is motivated by the fact that in such a 
situation the stress in the channel center is zero and 
linearly increases with distance from the middle of 
the channel. One thus expects a solid-like behavior 
in a central region (defined by a{z) < cr yie i d ; z de- 
noting the transverse coordinate measured from the 
channel center) and a fluid-like behavior in the re- 
gion between this central part and the channel walls 



(cr(z) > Oyieid). As shown in Fig. [9l this expecta- 
tion is nicely born out. An interesting consequence 
of this behavior is that, if for a given driving force 
(pressure gradient) the channel width is too small 
to allow the formation of stresses above cr yield , the 
whole system will behave like a solid body sticking 
to the walls. 
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